{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee1bc1b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import shutil\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import onekey_algo.custom.components as okcomp\n",
    "from onekey_algo import get_param_in_cwd\n",
    "\n",
    "plt.rcParams['figure.dpi'] = 300\n",
    "\n",
    "task = get_param_in_cwd('task_column') or 'label'\n",
    "bst_model = get_param_in_cwd('sel_model') or 'LR'\n",
    "model_names = ['Clinic Signature', 'Rad Signature']\n",
    "labels = [task]\n",
    "label_data_ = pd.read_csv('clinic_sel.csv')\n",
    "label_data_ = label_data_.dropna(axis=0)\n",
    "\n",
    "ids = label_data_['ID']\n",
    "print(label_data_.columns)\n",
    "label_data = label_data_[['ID', 'label']]\n",
    "label_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0453b415",
   "metadata": {},
   "source": [
    "# 汇总训练集"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4dc5642",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "Clinic_results = pd.merge(pd.read_csv(f'./results/Clinic_{bst_model}_train.csv', header=0), label_data, on='ID', how='inner')\n",
    "Rad_results = pd.merge(pd.read_csv(f'./results/Rad_{bst_model}_train.csv', header=0), label_data, on='ID', how='inner')\n",
    "\n",
    "ALL_results = pd.merge(Clinic_results, Rad_results, on='ID', how='inner')\n",
    "ALL_results.columns = ['ID', '-0', 'Clinic_Sig', task, '-00', 'Rad_Sig', '-l']\n",
    "\n",
    "ALL_results = ALL_results.dropna(axis=1)\n",
    "ALL_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c6b536f",
   "metadata": {},
   "source": [
    "# COX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e459cac",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines import CoxPHFitter\n",
    "\n",
    "cox_data = pd.merge(ALL_results[['ID', 'Clinic_Sig', 'Rad_Sig']], \n",
    "                    label_data_[['ID', 'label', 'duration']], on='ID', how='inner')\n",
    "\n",
    "cph = CoxPHFitter(penalizer=0.1)\n",
    "cph.fit(cox_data[[c for c in cox_data.columns if c != 'ID']], duration_col='duration', event_col='label')\n",
    "cph.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "061e0803",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components import nomogram\n",
    "\n",
    "cox_data.to_csv('results/nomo.csv', index=False)\n",
    "\n",
    "nomogram.nomogram(cox_data, duration='duration', result='label', columns=['Clinic_Sig', 'Rad_Sig'],\n",
    "                  survs=[6, 12*3, 12*5], surv_names=['half year survival','3 year survival','5 year survival'], with_r=False)\n",
    "\n",
    "# nomogram.risk_nomogram(cox_data, result='survival', columns=['Clinic_Sig', 'Rad_Sig'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7b03b49",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines.utils import concordance_index\n",
    "\n",
    "c_index_list = [[cph.concordance_index_, \n",
    "                 concordance_index(cox_data['duration'], -cox_data['Rad_Sig'], cox_data['label']),\n",
    "                 concordance_index(cox_data['duration'], -cox_data['Clinic_Sig'], cox_data['label']),\n",
    "                 'Train']]\n",
    "pd.DataFrame(c_index_list, columns=['Nomogram-Cox', 'Rad_Sig', 'Clinic_Sig', 'Cohort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "221f9f0a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pred_column = [f'{task}-0', f'{task}-1']\n",
    "gt = [np.array(d) for d in [Clinic_results[labels], Rad_results[labels]]]\n",
    "pred_train = [np.array(d) for d in [Clinic_results[pred_column], Rad_results[pred_column]]]\n",
    "okcomp.comp1.draw_roc(gt, pred_train, labels=model_names,  \n",
    "                      title=f\"Model AUC\")\n",
    "plt.savefig(f'img/train_auc.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a76a3eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.metrics import analysis_pred_binary\n",
    "metric = []\n",
    "for mname, y, score in zip(['Clinic Signature', 'Rad Signature'], gt, pred_train):\n",
    "    # 计算验证集指标\n",
    "    acc, auc, ci, tpr, tnr, ppv, npv, precision, recall, f1, thres = analysis_pred_binary(y, score)\n",
    "    ci = f\"{ci[0]:.4f} - {ci[1]:.4f}\"\n",
    "    metric.append((mname, acc, auc, ci, tpr, tnr, ppv, npv, precision, recall, f1, thres, f\"Shengyi\"))\n",
    "pd.DataFrame(metric, index=None, columns=['Signature', 'Accuracy', 'AUC', '95% CI', 'Sensitivity', 'Specificity', \n",
    "                                          'PPV', 'NPV', 'Precision', 'Recall', 'F1','Threshold', 'Cohort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5dd2842",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.delong import delong_roc_test\n",
    "\n",
    "results = pd.merge(Rad_results, Clinic_results, on='ID', how='inner')\n",
    "results.columns = ['ID', '-0', 'Rad', task, '-00', 'Clinic', '-l']\n",
    "\n",
    "delong = []\n",
    "delong.append([delong_roc_test(results[task], results[f'Rad'], results[f'Clinic'])[0][0], 'Train'])\n",
    "pd.DataFrame(delong, columns=['Rad Vs Clinic', 'cohort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9854a2a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.comp1 import plot_DCA\n",
    "plot_DCA([ALL_results[f'Clinic_Sig'], ALL_results[f'Rad_Sig']], \n",
    "         ALL_results[task], title=f'Model for DCA', labels=model_names)\n",
    "plt.savefig(f'img/train_dca.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac5fb028",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.comp1 import draw_calibration\n",
    "draw_calibration(pred_scores=pred_train, n_bins=10,\n",
    "                 y_test=gt, model_names=model_names)\n",
    "plt.savefig(f'img/train_cali.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89d1e859",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components import stats\n",
    "\n",
    "hosmer = []\n",
    "hosmer.append([stats.hosmer_lemeshow_test(y_true, y_pred[:,1], bins=15) \n",
    "              for fn, y_true, y_pred in zip(model_names, gt, pred_train)])\n",
    "pd.DataFrame(hosmer, columns=model_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a8fb4ee",
   "metadata": {},
   "source": [
    "# 测试集-TCGA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30ff5baa",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "Clinic_results = pd.merge(pd.read_csv(f'./results/Clinic_{bst_model}_test.csv', header=0), label_data, on='ID', how='inner')\n",
    "Rad_results = pd.merge(pd.read_csv(f'./results/Rad_{bst_model}_test.csv', header=0), label_data, on='ID', how='inner')\n",
    "\n",
    "ALL_results = pd.merge(Clinic_results, Rad_results, on='ID', how='inner')\n",
    "ALL_results.columns = ['ID', '-0', 'Clinic_Sig', task, '-00', 'Rad_Sig', '-l']\n",
    "\n",
    "ALL_results = ALL_results.dropna(axis=1)\n",
    "ALL_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06093016",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines import CoxPHFitter\n",
    "from lifelines.statistics import logrank_test\n",
    "from lifelines import KaplanMeierFitter\n",
    "\n",
    "cox_data = pd.merge(ALL_results[['ID', 'Clinic_Sig', 'Rad_Sig']], \n",
    "                    label_data_[['ID', 'label', 'duration']], on='ID', how='inner')\n",
    "\n",
    "c_index = cph.score(cox_data[[c for c in cox_data.columns if c != 'ID']], scoring_method=\"concordance_index\")\n",
    "# y_pred = cph.predict_median(cox_data[[c for c in cox_data.columns if c != 'ID']])\n",
    "# cox_data = pd.concat([cox_data, y_pred], axis=1)\n",
    "# mean = cox_data.describe()['duration']['mean']\n",
    "# cox_data['HR'] = cox_data[0.5] < mean\n",
    "\n",
    "y_pred = cph.predict_partial_hazard(cox_data[[c for c in cox_data.columns if c != 'ID']])\n",
    "cox_data = pd.concat([cox_data, y_pred], axis=1)\n",
    "cox_data['HR'] = cox_data[0] > 1\n",
    "\n",
    "dem = (cox_data[\"HR\"] == True)\n",
    "results = logrank_test(cox_data['duration'][dem], cox_data['duration'][~dem], \n",
    "                       event_observed_A=cox_data['label'][dem], event_observed_B=cox_data['label'][~dem])\n",
    "p_value = results.p_value\n",
    "kmf = KaplanMeierFitter()\n",
    "plt.title(f\"C-index:{c_index:.4f}, p_value={p_value}\")\n",
    "if sum(dem):\n",
    "    kmf.fit(cox_data['duration'][dem], event_observed=cox_data['label'][dem], label=\"High Rish\")\n",
    "    kmf.plot_survival_function(color='r')\n",
    "if sum(~dem):\n",
    "    kmf.fit(cox_data['duration'][~dem], event_observed=cox_data['label'][~dem], label=\"Low Risk\")\n",
    "    kmf.plot_survival_function(color='g')\n",
    "plt.savefig(f'img/km_TCGA.svg', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0f35b5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines.utils import concordance_index\n",
    "\n",
    "c_index_list.append([cph.concordance_index_, \n",
    "                     concordance_index(cox_data['duration'], -cox_data['Rad_Sig'], cox_data['label']),\n",
    "                     concordance_index(cox_data['duration'], -cox_data['Clinic_Sig'], cox_data['label']),\n",
    "                     'Test'])\n",
    "pd.DataFrame(c_index_list, columns=['Nomogram', 'Rad_Sig', 'Clinic_Sig', 'Cohort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b28d8126",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pred_column = [f'{task}-0', f'{task}-1']\n",
    "gt = [np.array(d) for d in [Clinic_results[labels], Rad_results[labels]]]\n",
    "pred_train = [np.array(d) for d in [Clinic_results[pred_column], Rad_results[pred_column]]]\n",
    "okcomp.comp1.draw_roc(gt, pred_train, labels=model_names,  \n",
    "                      title=f\"Model AUC\")\n",
    "plt.savefig(f'img/test_auc_TCGA.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b7a2896",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.metrics import analysis_pred_binary\n",
    "for mname, y, score in zip(['Clinic Signature', 'Rad Signature'], gt, pred_train):\n",
    "    # 计算验证集指标\n",
    "    acc, auc, ci, tpr, tnr, ppv, npv, precision, recall, f1, thres = analysis_pred_binary(y, score)\n",
    "    ci = f\"{ci[0]:.4f} - {ci[1]:.4f}\"\n",
    "    metric.append((mname, acc, auc, ci, tpr, tnr, ppv, npv, precision, recall, f1, thres, f\"TCGA\"))\n",
    "pd.DataFrame(metric, index=None, columns=['Signature', 'Accuracy', 'AUC', '95% CI', 'Sensitivity', 'Specificity', \n",
    "                                          'PPV', 'NPV', 'Precision', 'Recall', 'F1','Threshold', 'Cohort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "998ca8a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.delong import delong_roc_test\n",
    "\n",
    "results = pd.merge(Rad_results, Clinic_results, on='ID', how='inner')\n",
    "results.columns = ['ID', '-0', 'Rad', task, '-00', 'Clinic', '-l']\n",
    "\n",
    "delong.append([delong_roc_test(results[task], results[f'Rad'], results[f'Clinic'])[0][0], 'Test'])\n",
    "pd.DataFrame(delong, columns=['Rad Vs Clinic', 'cohort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1b1acef",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.comp1 import plot_DCA\n",
    "plot_DCA([ALL_results[f'Clinic_Sig'], ALL_results[f'Rad_Sig']], \n",
    "         ALL_results[task], title=f'Model for DCA', labels=model_names)\n",
    "plt.savefig(f'img/test_dca_TCGA.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62d6e04b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components.comp1 import draw_calibration\n",
    "draw_calibration(pred_scores=pred_train, n_bins=8,\n",
    "                 y_test=gt, model_names=model_names)\n",
    "plt.savefig(f'img/test_cali_TCGA.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6c8294d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from onekey_algo.custom.components import stats\n",
    "\n",
    "hosmer.append([stats.hosmer_lemeshow_test(y_true, y_pred[:,1], bins=8) \n",
    "              for fn, y_true, y_pred in zip(model_names, gt, pred_train)])\n",
    "pd.DataFrame(hosmer, columns=model_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e7796bb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
